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g Abstract 

Using integration by parts on Gaussian space we construct a Stein Unbiased 
Risk Estimator (SURE) for the drift of Gaussian processes using their local and 
CN occupation times. By almost-sure minimization of the SURE risk of shrinkage 

estimators we derive an estimation and de-noising procedure for an input signal 
perturbed by a continuous-time Gaussian noise. 
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1 Introduction 

Let X be a Gaussian random vector on M"^ with unknown mean m and known covari- 
ance matrix cr'^Id under a probability measure P^- 

It is well-known [13j that given (7 : ^ M"^ a sufficiently smooth function, the mean 
square risk ||X -|- g{X) — mW^^ of X -|- g{X) to m can be estimated unbiasedly by 

d d 

SURE :=a2ci + ^(?,(X)2 + 2^V,(7(X), (1.1) 

i=l i=l 
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from the identity 



[||X + g{X) - m\\l,] = a^d + E„ 



i=l 



1=1 



(1.2) 



which is obtained by Gaussian integration by parts under P^. The estimator (1.1) 
which is independent of m, is called the Stein Unbiased Risk Estimate (SURE). 



When {g^)x^K is a family of functions it makes sense to almost surely minimize the 



Stein Unbiased Risk Estimate (1.1) of (7 with respect to the parameter A. This point 



of view has been developed by Donoho and Johnstone [ij for the design of spatially 
adaptive estimators by shrinkage of wavelet coefficients of noisy data via 

X + g\X) = \r^{X/\), 

where 77 (x) is a threshold function. 



In this paper we construct a Stein type Unbiased Risk Estimator for the deterministic 
drift {ut)tm+ of a one dimensional Gaussian processes (Xt)tg[o,T] via an extension of 



the identity (1.2) introduced in |T0], [9] on the Wiener space. For example, given a{t) 
and \{t) two functions given in parametric form, the SURE risk of the estimator 

Xt + iT\X,) = a{t) + X{t)vs (^^i^) ' ^ e [0, T], 



where r]H is the hard threshold function ( 5.1[ ) below, is given by 



SURE(X + ^°'^(X)) = T + 



iXt-a{t)y 



-f{t,t) {\Xt-a{t)\<X./^)} 



dt + 2Xir - 2L 



T5 



where 7(s,t) = Cov{Xs, Xt), < s,t < T, denotes the covariance of {Xt)te[o,T] and 
i^, respectively denote the local and occupation time of 



(|Xi-a(t)|/v/^(M))tG[o. 



T], 



cf. Proposition 5.1 We apply this technique to de-noising and identification of the 
input signal in a Gaussian channel via the minimization of SURE {X + C,'^'^{X)). This 
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yields in particular an estimator of the drift of Xt from the estimation of a{t), and 
an optimal noise removal threshold from the estimation of A. This approach differs 
from classical signal detection techniques which usually rely on likelihood ratio tests, 
cf e.g. [S], Chapter VI. It also requires an a priori hypothesis on the parametric form 
of a{t). 

We proceed as follows. In Section [2] we recall our framework of functional estimation 
of drift trajectories. In Section [3] we derive Stein's unbiased risk estimate for the 
estimation of the drift of Gaussian processes. In Section |4] we discuss its application 
to soft thresholding for Gaussian processes using the local time and obtain an upper 
bound for the risk of such estimators. We also show the existence of an optimal 
parameter and the smoothness of the risk function. In Section |5] we consider the case 
of hard thresholding. In Section [6] we consider several numerical examples where a{t) 
is given in parametric form. In Section [7] we recall some elements of stochastic analysis 
of Gaussian processes. 

2 Functional drift estimation 

In this section we recall the setting of functional drift estimation to be used in this pa- 
per. Given T > we consider a real-valued centered Gaussian process X = {Xt)t£[o,T] 
with non-vanishing covariance function 

-f{s,t)=-E[X,Xt], s,te[0,T], 

on a probability space {Q, JF, P), where {J-')te[o,T] is the filtration generated by {Xt)te[o,T] 
Assume that under a probability measure we observe the paths of {Xt)t,= [o^T] de- 
composed as 

Xt = ut + x^, te[o,T], 

where u = {ut)te[o,T] is a square integrable jF-adapted process and (X")tg[o,T] is a 
centered Gaussian process with covariance 

7(s,^)=E„[X,^X«], 0<s,t<T, 
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(2.1) 



where denotes the expectation under P^. Given a continuous time observation of 
the process {Xt)t<^[o^T] we will propose estimators of the unknown drift function u. 

Definition 2.1. The risk of an estimator C, '■= {C,t)te[o,T] to u is defined as 

Jo 

where fi is a positive measure on [0,T]. 

Examples of risk measures /i include the Lebesgue measure and 

n 

(x{dt) = ^aiSt^{dt), ai,...,an>0, 

i=l 

in which case the risk of the estimator is computed from discrete values of the sample 
path observed at times ti, . . . , t„, n > 1. 

Definition 2.2. A drift: estimator {C,t)t£[o,T] is called unbiased if 

E46]=E4Mt], te[o,T], 

for all square-integrable J^t-o^dapted process {ut)te[Q,T], where {J^t)t(^[o.T] is the filtration 
generated by {Xt)t(i[o,T\- 

In the sequel we will consider the canonical process {Xt)t(z[o,T] as an unbiased estimator 
u := (Xt)ig[o,T] of its own drift (ut)te[o,T] under P^j, with risk 



R(7,/i,M) := E„ 



\Xt - ut\^n{dt) 



■y{t, t)fi{dt) 



Recall that the estimator u = (Xt)tg[o,T] is minimax i.e. 



R(7, /i, u) = inf sup E^, 

5 ven 



16 - Vtl'^Kdt) 



cf. Proposition 3.2 of [lO]. In addition, when {Xt)t£[o,T] has independent increments 
and (Mt)tg[o,T] £ L'^{^ x [0,T],P„ (g) /i) is square-integrable and adapted, then for any 
adapted and unbiased estimator ^ the Cramer-Rao bound 



E„ 



/ - utl"^ ij{dt) 
Jo 



> R(7,/i,M) 



(2.2) 



holds for any unbiased and adapted estimator {C,)te[o,t] of {u(t))te[o,T] G L^{Q x 
[0, T], P„ (g) fi) and is attained by u, cf. Proposition 4.3 of [lOj, hence u = {Xt)t<=io,T] is 
an efficient estimator of its own drift u. 
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3 Stein's unbiased risk estimate 



Instead of using the minimax estimator u we will estimate the drift of (Xj^gjo^rj by the 
almost sure minimization of a Stein Unbiased Risk Estimator for Gaussian processes, 



constructed in the next proposition by analogy with (1.1). In the next proposition 
we use the gradient operator Dt whose definition and properties are recalled in the 



appendix, cf. Definition 7.2 and Lemma 7.3 



Proposition 3.1. For any (6)t6[o,r] ^ L'^{^ x [0, T], P„ such that ^ Dom(V), 
t G [0,T], and (A6)te[o,T] G L^(fi x [0,T],P„ (g) /i), the quantity 



SURE^(X + e) := R(7,/^,«) + l|e||i.([o,T] 



d^l) 



Dtitii{dt) 



is an unbiased estimator of the mean square risk ||X + ^ — m"^ 



L2([0,T],dAt)- 



Proof. From Lemma [7. 3| we have 



XtfKdt) 
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[/ 

Jo 


x: + 6 


fi{dt) 



lief 



R(7,/i,M) + E„ llell 
E„ [SURE^(X + OI- 



LS ([o,T],dA») 
+ 2E, 



+ 2E, 

T 

DtiHdt) 







(3.1) 



□ 



Unlike the pointwise mean square risk ||X + ^ — m||^2([ot] d/^)' SURE risk estimator 
does not depend on the estimated parameter u. 

Given a family (^'*')AeA of estimators indexed by a parameter space A, we consider the 
estimator X + that almost-surely minimizes the SURE risk, with 

A* = argmin;,g^SURE^(X + ^^). 

For all values of A the SURE risk of the estimator X + C,^* improves on the mean 
square risk of X + 
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Precisely for all u G A we have 

E„[SURE^(X + e^*)] < E„[SURE^(X + r)] 



\r-u\\u 



inf E„ 

A 



11^ ~ '"llL2([0,T],/i) 



In the sequel we will apply the above to a process (6)t6[o,T] given as a funtion = 
C,t{Xt) of Xt, t G [0,T]. In particular we will discuss estimation and thresholding for 
estimators of the form 

Xt + = a{t) + X{t)r^ (^^A(f^) ' ^^'^^ 

where : M — > M is a threshold function with support in (—00, — 1] U [1, 00). 

In particular we will apply our method to the joint estimation of parameters a, A, 



successively in case a{t) = a, a{t) = at, and X{t) = Aa/tIvO- 

4 Soft threshold 

In this section we construct an example of SURE shrinkage by soft thresholding in 



the framework of Proposition 3.1, with application to identification and de-noising in 



a Gaussian signal. In case rj is the soft threshold function 

risiy) = sign{y){\y\-iy, y eR, (4.1) 



the function ^"'^ in (3.2) becomes 



at 



(x) = — sign(a; — a{t)) min(A(t), \x — a{t)\), x G M, 
where A(t) > is a given level function. 
Proposition 4.1. We have F-a.s 

SVRE ^{X + C^{X)) (4.2) 

= R(7,/i,u)+ /" \Xt-a{t)\^ AX\t)ii{dt)-2 [ l||x,-a(t)|<A(t)}7(i, i)/^(c?i)- 
Jo Jo 
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Proof. Since 



L{|x-a(t)|<A(i)}, 



we have 



T 



{t)\<\{t)}DtXtfi{dt) 

l{|Xt-a{i)|<A{i)}7(^, t)fl{dt), 



hence the conclusion from Proposition 3.1 



□ 



The risk associated to discrete observations (X^^, . . . , Xt„) can be computed via Propo- 
sition 



4.1 by choosing the risk measure (2.1), in which case Relation (4.2) becomes 



SURE {X + C\X)) 

n n 

= R{-f,fi,u) + ^\Xt^ -a(t,)r A A2(ti) -2 J^7(ti,ti)l{|x,^_„(t^)l<A(t,)}. 



i=l i=l 

which is analog to the finite dimensional SURE risk 

d 



SURE (X + g\X)) = d + J2i\^i\^ - 2#{^; l^^il < A} 



(4.3) 



of [3] . In the simulations of Section |6] we effectively use such risk measures when 
discretizing the signal. More precisely, when fj,{dt) = f{t)dt has a density f{t) with 
respect to the Lebesgue measure and 

n-l 



i=l 



Relation (4.2) shows that SURE^^(X + (X)) becomes a consistent estimator of 



the risk SURE^(X + (^)) as n goes to infinity. 



Taking 

H{dt) =-f'^{t,t)dt and X{t) = X^/-f{t,t), A > 0, tG[0,T], 
and letting 

hx.-.it)l<xV^)}dt (4.4) 
7 



denote the occupation time of the process 



up to time T in the set [—A, A], Proposition 4.1 yields the identity 

SURE^(X + r'^(X)) =r+ / A A)^rft-2L^. (4.5) 

As a consequence we obtain the following bound for the risk of the thresholding 
estimator X + 

Proposition 4.2. Assume that u G L^([0, T], (i/i) is a deterministic function and let 
fi{dt) := 'j{t,t)~^dt. Then for all fixed X > we have 

E„[||X+r''(X)-tz||i.([o,r],,^)] < (l+A^) (TA£\u{t)-a{t)\'fxidt)^+Til+\)e-'^. 
Proof. We have 

pT pT 



SURE f,{X + C'\X)) = T + {\Z^'^\AXfdt-2j^ r^.dt 



hence 

1E„[SURE^(X + r'\X))] < T(l + A^), 

and 

E4SURE^(X + e'\X))] < f 1 + -^l^] A A^ - 2P„(|Zf < \)dt 

Jo 



< (l + A^) / |M(t)-a(t)|V(c^t) + r(l + A2)e-V 
where we recall that from [3], Appendix 1, we have for every t in [0,T] that 



1 + E„[|Z°'^|^] A A^ - 2P„(|Zf'^| < A) < (1 + A^) e 



l{t,t) 



and we conclude from Proposition 3.1 □ 
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From this proposition it follows that SURE^(X + ^"''^(X)) is independent of large 
values \\u — a;||i2Qo^T']), while its growth at most as 1 + in A > 0. 



Since A i— SURE^(X + ^"''^(X)) in (4.5) is lower bounded by — T and equal to when 
A = 0, the optimal threshold 

A* := argmin;,SURE ^(X + C'^{X)) (4.6) 

exists almost surely in [0, oo). 

In addition we have the following proposition which important for the numerical search 
of an optimal parameter value. 

Proposition 4.3. The function A ^ SURE^(X + ^°''^(X)) is continuously differen- 
tiable. 

Proof. Letting 

A{s,t) = Var„(Zr^ - Zf'^) = 2 - 2 0<s,t<T, 

V7(s,s)7(t,t) 



under Condition (7.4), the local time 



of (|Z"''^|)fg[o,T] exists almost surely, cf. Section [7| and we have 

|^SURE,(X + r'\X)) = ^l\\Zr\AXfdt-2i'^ 



2^ 1^ ^{\x,-ait)\>x^)}dt - 



= 2A(T - L^) - 2£^, 

which is a continuous function of A since the covariance 7(s,t) does not vanish, cf. 
e.g. Theorem 26.1 of [H]. □ 
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-SURE^(X + r''W)|A=o = -21?., 



Consequently we have 

d_ 

hence A* > a.s. when is a.s. positive, which is the case for example when Xf is 
a Brownian motion, see Corollary 2.2 of page 240 of [12], Chapter VI. 



In practice we will compute A* numerically by minimization of A i— SURE^(X + 
^"''^(X)) over A in a range A = [0, C{T)] where C(T) is such that 

lim P„ ( sup \Z^'^\ < C{T)] = 1. 



T->oo 



te[o,T] 



This condition is analog to Condition (31) in [3] and allows us to restrict the range of 
A when searching for an optimal threshold. 

The function a{t) can be given in parametric form, in which case the parameters will 
be used to minimize SURE^(X + ^"'^{X)), cf. Section g 



5 Hard threshold 

Here we use the threshold function 

Vniy) = l/l{|y|>i}, 2/ e (5.1) 

hence 

-{\x-a(t)\<\y/^)r 



where A > is a level parameter. 



In finite dimensions [3] the SURE estimator (1.1) can not be computed due to the 
non-differentiability of rjH, however a deterministic optimal threshold equal to ■\/2Togci 
can be obtained by other methods, cf. Theorem 4 of [3]. 

In continuous time the situation is different due to the smoothing effect of the integral 
over time. In the next proposition we compute the SURE risk using the local time of 
Gaussian processes when fi{dt) = 'y~^{t,t)dt. 
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Proposition 5.1. We have F-a.s 

7(t,t) ^{\Xt-a{t)\<\y/^)f 



(5.2) 







Proof. Let (p G C^([— 1,1]), > be symmetric around the origin, such that 
J^^ (f){x)dx = 1, and let 

Let 

/oo 

denote the convolution of ^e^/:^^^) ^^^^ ^t'^^ with 



From the occupation time density formula (7.5) we have 

T 



T roo 




A / ^,{a-\)i^da 

(y)^{\.-y-ait)\<X^)}dydt^ 



J~oc 
oo 



oo 

T roo 




J-oc 
oo 



oo 

T roo 




'0 J-oo 

11 



which converges in L'^{Q,Fu) to 

- ^{\X,-a{t)\<X,/^)/^ 

as e tends to zero. □ 



6 Numerical examples 

In this section we assume that is a centered stationary Ornstein-Uhlenbeck process 
solution of 

dX^ = -aX^dt + adBt, t e [0, T], 

2 

with X^ ~ (o, and covariance function 'y{s,t) = ^e""'*"''', s,t G [0,T], for 
(T,a > 0. As a consequence of the following proposition we can take A = [0, \/2 logT] 
as parameter range when T is large. 

Proposition 6.1. Assume that ||«||l=°([o,oo)) < oo and ||m||l°°([o,oo)) < oo. Then for 
any r > 1 we have 



lim P„ sup \Zt\ < ^/2r\ogT | = 1. 
\te[o,T] 

Proof. From Theorem 1.1 of [14j (see also [7], Theorem 2.1 of [IT], and page 488) 
there exists a universal constants ci, C2 > such that for all A,T > 0, 

P„ I sup |^t|>Al < ciM(2aT,C2/A)^(A), 

\t&[0,T] 



where = J!^ e ^ ^'^dy/\/2TT and M(2aT, C2/A) is the maximal cardinal of all 

sequences S in [0, 2aT] such that 



I 7 7 II / ^ - ^ C2 _ 

I^t - ^s||L2(n) = a"' S,tE<b. 

Setting A = V2rlogT, r > 0, T > 1, and using the bound ^(A) < e'^^^^ /{XV2Tr) 
this yields, for all T large enough: 

P I sup \Zt\< ^2r logT I > 1 - c^T^-^ 
Ve[o,T] y VO 

which tends to 1 as T ^ 00 provided r > 1. □ 
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In the next figures we present some numerical simulations when the signal {Xt)te[o,T] 
is a deterministic function (u(t))fg[o,T] perturbed by a centered Ornstein-Uhlenbeck 
process, with parameters a = 0.5, a = 0.05, T = 1. 



We represent simulated samples path with the optimal thresholds obtained by soft 
thresholding, the de-noised signal after hard thresholding, and the corresponding risk 
function (a. A) ^ SURE^(X + ^"'^(X)) whose minimum gives the optimal param- 
eter value(s). The hard threshold function has not been used for estimation due to 



increased numerical instabilities linked to the simulation of the local time in (5.2). 



Simple thresholding 

Here we take Ut = 0.2 x max(0, sin(37rt)), X{t) 
signal around the level a{t) = 0, t E [0,T]. 



and we aim at de- noising the 




0.04 0.06 
threshold value 



Figure 6.1: Risk function A ^ SURE/,(X + 



From Figure 6.1 we estimate the optimal threshold to A*w^ = 0.018, after numerical 



minimization on a grid, which leads to the thresholding described in Figure 6.2 below. 
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Figure 6.2: Process trajectory 



Estimated trajectory 



Level detection and thresholding 

We apply our method to the joint estimation of parameters a, A, in case Ut = 0.3 + 
0.2sign(sin(27rt)) xmax(0, sin(37rt)), a(t) = a and X(t) = X^/^, i.e. we aim at detecting 
simultaneously the level a = 0.3 and the threshold X^y^ at which the noise can be 
removed. For this we have the following proposition that completes Proposition 4.3. 

□ 

Proposition 6.2. The function (a, A) i-^ SURE ^(X + is continuously dif- 

ferentiable. 



Proof. We have 
d 



da 



where i^' denotes the local time at level a of the process 




' ^ — 0.8 

0.4 threshold value 



Figure 6.3: Risk function (a, A) ^ SURE/,(X + (^)) 
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From Figure 6.3 we estimate the optimal threshold and shift parameters at A*^y7 = 
0.017 and a* = 0.30, which leads to the thresholding described in Figure 6.4 below. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0. 




0.2 0.3 0.4 0.5 0.6 0.7 0.8 



Figure 6.4: Process trajectory 



Estimated trajectory 



Figure 6.3 also shows that the values 0.5 and 0.1 are other candidates to an estimation 
of a. iQe values correspond to the extrema in the sample trajectory of Figure 6.4. 

□ 

Drift detection and thresholding 

We apply our method to the joint estimation of parameters a, A, in case Ut = 0.3t + 
0.2sign(sin(27rt)) x max(0, sin(37r)f:)), a(t) = at, and X(t) = X^/j, i.e. we aim at 
locating noise with threshold around a line of slope a = 0.3. Analogously to 
Propositions 4.3 and 6.2 we have the following result. 
□ □ 

Proposition 6.3. The function (a, A) 1-^ SURE ^(X + is continuously dif- 

ferentiable. 

Proof. We have 

d 



SURE^(X + r'^(X)) = -2 



da 



Xf-at 



where denotes the local time at level a of the process 
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Figure 6.5: Risk function (a, A) ^ SURE^(X + 

The optimal threshold and slope parameters are numerically estimated at A*^/7 = 
0.0093 and a* = 0.294. 




Figure 6.6: Process trajectory Estimated trajectory 



The threshold and slope and actually slightly underestimated, as the larger noise at 
the right end of the slope line has been interpreted as being part of the signal. 

7 Appendix 

In this section we review three aspects of stochastic analysis for Gaussian processes, 
including local time and the Malliavin calculus calculus. 
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Malliavin calculus on Gaussian space 

Here we recall some elements of the Malliavin calculus on Gaussian space for the 
centered Gaussian process {Xt)t(^[o,T], see e.g. [6]. Let be a finite Borel measure on 
[0, T] and let T the operator defined as 

iTg)it)= [ gisMs,t)fx{ds), t e [0,T], 
Jo 

on the Hilbert space H of functions on [0,T] with the inner product 

{h,g)H = {h,Tg)L2(^[Q^T],di^)- 
The process {Xt)te[o,T] can be used to construct an isometry X : H ^ L'^{Q, JF, P) as 

X{h) = [ X,h{s)fx{ds), heH. 
Jo 

Then {X{h) : h G H} is an isonormal Gaussian process on H, i.e. a family of 
centered Gaussian random variables satisfying 

nX{h)X{g)] = {h,g)H, Kg EH. 

For any orthonormal basis {hk)km of L^([0, T], d/i), we have the Karhunen-Loeve 
expansion 

oo 

Xt = Y,hk{t)X{hk), te[0,T]. (7.1) 

A:=0 

Let now S denote the space of cylindrical functionals of the form 

F = /„(X"(/ii),...,X"(/i„)), (7.2) 

where /„ is in the space of infinitely differentiable rapidly decreasing functions on M", 
n > 1. 

Definition 7.1. The H -valued Malliavin derivative is defined as 

n 

VtF = J2 Ht)^^fn (X"(/li), . . . , X\K)) , 
i=l 



for F E S of the form (7.2). 
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It is known that V is closable, cf. Proposition 1.2.1 of [H], and its closed domain will 
be denoted by Dom(V). 

Definition 7.2. Let Dt he defined on F & Dom (V) as 

DtF:={TVF){t), tG [0,T]. 

Let 6 : Lf^{Q] H) L^(r2,P„) denote the closable adjoint of V, i.e. the divergence 
operator under P„, which satisfies the integration by parts formula 

]E^[F6{v)] = ]E^[{v,VF)h], FeDom(V), veBom{6), (7.3) 

where 1E„ denotes the expectation under P^, with the relation 

S{hF) = FX{h)-{h,VF)H, 

cf. [6], for F G Dom ( V) and h ^ H such that hF G Dom (5). The next lemma will 



be needed in Proposition |3.1| below to establish Stein's Unbiased Risk Estimate for 
Gaussian processes. 

Lemma 7.3. For any F G Dom (V) and u E H we have 

E„[FXr] =E„[AF], t G [0,T]. 

Proof. We have 

oo 

-E^IFX^] = ^/i,,(t)E„[FX"(/ifc)] 

fc=0 

oo 

fc=0 

oo 

= 5^/ifc(t)E„[(/ifc,VF)^,] 



fc=0 

oo 



Y,hk{t) E„[(/ifc,rVF)z.2([o,T],M)] 

fc=0 

E„ [(F VF) (t)] , F G Dom ( V) , t G [0, T] . 

□ 
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Note that since u & H we have VsXt{h) = VsX"(/i) = h{s) and 
DtXt = {TVXt){t) 







oo „T 
l,_n "'0 



fc=0 

oo 



k=0 

= ^{t,t), te[o,T]. 

Local time of Gaussian processes 

Given {Zt)te[o,T] a Gaussian process let 

A{s,t) =VaT{Zt- Zs), 0<s,t<T, 

and denote by 

:= / l!Zt<ndt 



the occupation time of (Zt)tg[o,T] up to T in the set (— cxd,A]. 

Recall that a classical result of Berman see Theorem 21.9 of [5], shows that if 

[ [ A-\s, t)dsdt < oo, (7.4) 
Jo Jo 

then for any A G M the local time 

f\ ._ _^rA 
- QX^T 

of (Zj)tg[o,T] at the level A exists and the occupation time density formula 

f f{Zt)dt= [ f{\)i'^d\ (7.5) 
Jo Jr 

holds for every positive measurable function / on M. The local time of \Zt\ is given 
by = i^"" + and the related occupation time formula can be obtained under the 
same condition from the relation 

/•T /'OO /"OO 

I fi\Zt\)dt= / f{\a\)tTda= / fia)tTda. 

J-oo Jo 
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